Tema 10 - Modelización con tidymodels

Técnicas para ‘Big Data’ en Economía - Curso 2023/24
Universidad de Alicante

Pedro Albarrán

Dpto. de Fundamentos del Análisis Económico. Universidad de Alicante

Alberto Pérez Bernabeu

Proceso de modelización

  • tidymodels es una colección de paquetes para el proceso de modelización (NO implementa modelos) con los principios de tidyverse
install.packages("tidymodels")
library(tidymodels)
  • Otros paquetes “similares”: mlr3, caret, H2O

Pre-procesado: partición inicial

  • dplyr transforma los datos para adecuarlos a la modelización, pero tidymodels permite transformaciones específicas para un modelo concreto

  • initial_split(): particionar los datos en prueba y entrenamiento.

library(mosaicData)
set.seed(9753)
RailTrail_part <- RailTrail %>% initial_split(prop = .8)

censo <- read_csv("data/census.csv") %>% mutate(income = factor(income))
set.seed(7482)
censo_part <- censo %>% initial_split(prop = .8)
  • Las funciones training() y testing() acceden a cada submuestra
RailTrail_entren <- RailTrail_part %>% training()
RailTrail_prueb  <- RailTrail_part %>% testing()

Pre-procesado: recetas

  • Las recetas definen las transformaciones a aplicar

  • recipe() tiene como primer argumento una fórmula (rol de las variables)

  • Luego se añaden pasos con step_

    • step_filter(), step_arrange(), step_rm(), etc.

    • step_naomit()

    • step_impute_mean(), step_impute_linear(), step_impute_knn()

  • Se aplican a una variable, todas o un subconjunto con all_outcomes(), all_predictors(), all_numeric(), all_nominal(), contains(), etc.

    • step_dummy(all_predictors(), -all_numeric())

Pre-procesado: recetas (cont.)

  • Se pueden añadir términos polinomiales con step_poly(), interacciones de variables con step_interact(), discretizar variables con step_cut(),etc.

  • Se centran variables con step_center() o se estandarizan con step_scale()

  • Se crea un objeto de receta (y luego combinarlo con otras partes del proceso)

receta_lm1 <- training(RailTrail_part) %>%            # Datos: NO crucial
  recipe(volume ~ cloudcover + precip + avgtemp) %>%
  step_poly(avgtemp, degree = 6) %>%                  
  step_center(all_predictors(), -all_nominal()) %>%
  step_scale(all_numeric(), -all_outcomes()) 

receta_logit1 <- training(censo_part) %>%
  recipe(income ~ age + education + race + sex + capital_gain)
  • También se podría preparar la receta con prep() y aplicar las transformaciones con juice() o bake()

Definir el modelo

  • tidymodels define un modelo con una interfaz unificada para distintas bibliotecas (con otros nombres de argumentos)
modelo_lm1     <- linear_reg(mode= "regression", penalty = 0) %>%
                    set_engine("lm") 

modelo_glmnet1 <- linear_reg(penalty = 0)  %>% set_mode("regression") %>% 
                        set_engine("glmnet")
modelo_logit1 <- logistic_reg(mode= "classification", penalty = 0) %>% set_engine("glm")
  • Se podría estimar (entrenar) el modelo con fit()

Flujos de trabajo: workflow()

  • Combina preprocesado y definición del modelo en un único objeto de flujos
lm1_flujo <- workflow() %>%
  add_recipe(receta_lm1) %>%      
  add_model(modelo_lm1)
  • Un flujo existente se modifica con update_recipe() , update_model(), etc.
  • Se prepara todo en una única llamada de fit()
lm1_flujo_est <- lm1_flujo %>% fit(data = RailTrail_part %>% training()) 

logit1_flujo_est <-  workflow() %>% add_recipe(receta_logit1) %>%      
                        add_model(modelo_logit1) %>% fit(censo_part %>% training())
  • Este objeto contiene tanto la receta para transformar los datos como el modelo estimado para mostrar los resultados o predecir

Flujos de trabajo: workflow() (cont.)

  • Se puede extraer la receta para aplicar la transformación a unos datos (p.e., comprobamos que tienen varianza 1 con var())
lm1_flujo_est %>% extract_recipe() %>% bake(RailTrail_part %>% training()) 
lm1_flujo_est %>% extract_recipe() %>% bake(RailTrail_part %>% testing())
  • También se pueden extraer los resultados de la estimación. Con funciones de broom se convierten a tibbles para trabajar con ellos (p.e., kable() en un documento)
lm1_flujo_est %>% extract_fit_parsnip() %>% tidy()      # resultados de la estimación
lm1_flujo_est %>% extract_fit_parsnip() %>% glance()    # otros detalles de la estimación
  • broom::augment() calcula predicciones del modelo, residuos, etc.
modelo <- lm1_flujo_est %>% extract_fit_parsnip()
modelo$fit %>% augment()      

Predicción

  • La función predict() (de parsnip) predice valores numéricos, clases, probabilidad de cada categoría, intervalos de confianza, etc.

  • Devuelve un tibble que podemos añadir a los datos con bind_cols()

 lm1_flujo_est %>% 
  predict(new_data = RailTrail_prueb) %>% 
  bind_cols(RailTrail_prueb %>% select(volume))     # Variable predicha .pred

logit1_flujo_est %>% 
  predict(censo_part %>% testing())                # clase predicha .pred_class

logit1_flujo_est %>% 
  predict(censo_part %>% testing(), type = "prob")    # probabilidades de cada categoría
  • Notad que se procesan automáticamente los datos donde se predice

Validación del Modelo

  • Dado un tibble con valores/clases observados (truth) y predichos (estimate), se calcula una métrica (rmse(), rsq(), accuracy(), etc.) o varias (metrics())
lm1_flujo_est %>% predict(RailTrail_prueb) %>% 
  bind_cols(RailTrail_prueb) %>%  
  metrics(truth=volume, estimate= .pred)           #   mae(truth=volume, estimate= .pred)

logit1_flujo_est %>%  predict(censo_part %>% testing()) %>% 
   bind_cols(censo_part %>% testing()) %>%   
  conf_mat(truth=income, estimate= .pred_class)           # clase predicha=
                                                          #  mayor probabilidad predicha

mis_metricas <- metric_set(accuracy, spec, sens)             # métricas específicas
logit1_flujo_est %>%  predict(censo_part %>% testing()) %>% 
  bind_cols(censo_part %>% testing()) %>%
  mis_metricas(truth=income, estimate= .pred_class)   

Validación del Modelo (cont.)

  • Si predecimos probabilidades, se pueden obtener la curva ROC y la AUC
logit1_probs <- logit1_flujo_est %>% 
                    predict(censo_part %>% testing(), type = "prob") %>%
                    bind_cols(censo_part %>% testing()) 

logit1_probs %>% roc_auc(income, `.pred_<=50K`) 

logit1_probs %>% roc_curve(income, `.pred_<=50K`) %>%  autoplot()
  • Con más de dos clases, se predice la probabilidad de cada clase y la clase predicha es la más frecuente

    • La matriz de confusión es similar, pero de dimensiones \(\small k \times k\)

    • La accuracy sigue teniendo la misma interpretación

    • La ROC-AUC se calcula para cada clase frente a las demás

El proceso de validación cruzada

  • vfold_cv() crea las particiones en los datos sin procesar
set.seed(9753)
RailTrail_cv <- RailTrail %>% vfold_cv(v=10) # objeto de 10 grupos
  • Se puede acceder a los datos de entrenamento y prueba de cada uno de ellos con analysis() y assessment(), respectivamente
RailTrail_cv$splits[[1]] %>% analysis() %>% dim()
RailTrail_cv$splits[[1]] %>% assessment() %>% dim()
  • fit_resamples(), similar a fit(), sobre un flujo de trabajo ya definido y el objeto completo de remuestras de valicación cruzada …
lm1_flujo_cv_est <- lm1_flujo  %>% 
                        fit_resamples(RailTrail_cv)

El proceso de validación cruzada (cont.)

  • … y el objeto creado contiene los valores de las métricas
lm1_flujo_cv_est %>% collect_metrics()      # promedio sobre 10 iteraciones
lm1_flujo_cv_est$.metrics %>% bind_rows()   # valores en cada iteración
  • Se pueden cambiar varias opciones del proceso
lm1_flujo_cv_est <- lm1_flujo  %>% 
                        fit_resamples(
                          resamples = RailTrail_cv, 
                          metrics   = metric_set(rmse, mae),
                          control   = control_resamples(save_pred = TRUE)
                        )

Selección de hiperparámetros: tuning

  • Algunos parámetros, denominados hiperparámetros, no pueden aprenderse directamente durante el entrenamiento del modelo (ej., \(\small \lambda\) en LASSO)
  • Proceso de ajuste (tuning):

    • en una parte de la muestra de entrenamiento, se estiman los parámetros dado un valor del hiper-parámetro

    • en la otra parte, se mide el error asociado a ese hiper-parámetro para validar el mejor valor

    • se elige el valor con mejor métrica de error en validación

\(\Downarrow\)

\(\hspace{0.1cm}\)

Proceso de tuning

  • Se pueden identificar los hiperparámetros a ajustar en la receta y/o el modelo
modelo_LASSO <- linear_reg(mode= "regression", 
                           penalty = tune() ) %>%
                    set_engine("glmnet") 

flujo_LASSO_tuning <- workflow() %>%
  add_recipe(receta_lm1) %>% 
  add_model(modelo_LASSO)

flujo_LASSO_tuning %>% parameters()
  • Se establecen las combinaciones de valores sobre las que se buscará con grid_random(), grid_max_entropy() o grid_regular()
LASSO_grid <- grid_regular(penalty(range = c(0, 15), trans = NULL), levels = 51)
                                   # rango,                 # número de valores

Proceso de tuning (cont.)

  • Se usa tune_grid() de forma a similar a fit_resamples() en las remuestras por Validación Cruzada de la muestra de entrenamiento
set.seed(9753)
RailTrail_entren_cv <- RailTrail_entren %>% vfold_cv(v=10)

LASSO_tuned <- flujo_LASSO_tuning %>% tune_grid( resamples = RailTrail_entren_cv, 
                                                 metrics   = metric_set(rmse, mae),
                                                 grid      = LASSO_grid               )
  • Podemos explorar visualmente el ajuste para distintos valores
LASSO_tuned %>% autoplot()

penalty <- LASSO_tuned %>% collect_metrics() %>% filter(.metric == "mae")
penalty %>% ggplot(aes(x=penalty, y=mean)) +  scale_x_log10() +
              geom_line() + geom_point(color="red") + 
              geom_errorbar(aes(ymin=mean-std_err, ymax=mean+std_err), color="gray")

Estimación del modelo completo

  • NOTA: debemos probar manualmente varios rangos y valores buscando la U

  • Se puede ver numéricamente el mejor o los cinco mejores candidatos: recordar que, por la variabilidad, varios valores son igualmente aceptables

LASSO_tuned %>% show_best("mae")
mejor_lambda <- LASSO_tuned %>% select_best("rmse")
  • Actualizamos el flujo de trabajo con un valor, ej., de select_best()
flujo_final <- flujo_LASSO_tuning %>% 
        finalize_workflow(mejor_lambda)  # finalize_workflow(parameters = list(penalty=8.5))
  • Debemos estimar el modelo en los datos completos de entrenamiento
flujo_final %>%  fit(data = RailTrail_entren) %>%  broom::tidy()

Finalizando y evaluando el modelo

  • Finalmente, podemos usar last_fit(): ajusta al modelo finalizado en los datos de entrenamiento y lo evalúa en los de prueba.
final_fit <- flujo_final %>% last_fit(RailTrail_part)

final_fit %>% collect_metrics()
final_fit %>% collect_predictions()


final_fit <- flujo_final %>% last_fit(split   = RailTrail_part,
                                      metrics = metric_set(rmse, mae))
final_fit %>% collect_metrics()